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1 .  Introduction 

In  recent  years,  there  has  been  strong  interest  in  applying  pseudospectral 
methods  to  various  flow  problems.  Some  examples  of  these  applications  are 
collected  in  a  recently  published  book  by  Voight,  Gottlieb,  and  Hussaini 
(1984).  One  of  the  current  areas  of  interest  is  the  calculation  of  compres¬ 
sible  flows  with  shock  waves.  Gottlieb,  Lustman,  and  Orszag  (1981)  have 
investigated  the  one-dimensional  shock  tube  problem  using  a  pseudospectral 
method  and  reported  the  ability  of  capturing  the  shock  within  one  grid  point. 
Gottlieb,  Lustman,  and  Streect  (1982)  have  reported  work  on  two  problems  in 
this  area.  The  first  is  the  solution  of  Euler  equations  for  oblique  shock 
reflection  from  a  flat  plate.  By  using  a  sparse  8x8  computational  grid,  they 
have  shown  that  it  is  possible  to  capture  the  shock  wave  within  one  grid  point 
and  that  the  treatment  of  boundary  conditions  is  extremely  crucial  to  the  con¬ 
struction  of  a  stable  scheme.  The  second  problem  is  the  solution  of  a  full 
potential  equation  for  transonic  flow  past  an  airfoil.  An  airfoil  is  mapped 
to  a  circle.  In  the  circumferential  direction,  a  Fourier  series  is  used,  while 
Chebyshev  expansion  is  used  in  the  radial  direction.  In-order  to  stabilize 
the  calculation,  an  artificial  viscosity  is  used  in  the  governing  equation. 
Gottlieb  et  al.  (1982)  has  shown  that  in  a  subsonic  flow,  a  highly  accurate 
solution  can  be  obtained  by  using  a  sparse  grid  of  32x8  grid  points.  In  the 
transonic  case  with  shocks,  the  shock  wave  spans  across  three  grid  points  and, 
hence,  is  not  accurate  enough  if  a  sparse  grid  is  used.  The  problem  chosen  by 
Gottlieb  et  al.  (1982)  for  the  solution  of  the  Euler  equations  is,  unfor¬ 
tunately,  not  a  good  one.  The  region  between  the  oblique  shock  is  a  constant 
state.  Aside  from  capturing  the  shock  wave,  their  research  yields  no  informa¬ 
tion  about  the  accuracy  of  the  method  in  the  region  of  smooth  variation  away 
from  the  shock  wave.  Their  work  on  potential  flows  has  shown  that  the 
incorporation  of  conventional  artificial  viscosity  into  the  spectral  method 
will  cause  deterioration  of  the  accuracy  and  thus  defeat  the  purpose  of  using 
a  spectral  method. 

In  the  present  work,  we  have  chosen  the  realistic  problem  of  transonic 
flow  over  an  airfoil  to  study  the  application  of  a  spectral  method  to  compres¬ 
sible  flows  with  shock  waves.  Part  of  the  findings  of  the  present  work  have 
been  reported  in  a  conference  paper  (Jou,  Jameson,  and  Metcalfe,  1983),  which 
is  included  in  this  report  as  an  appendix. 
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2 .  Basic  Approach 
2.1  Governing  Equations 

The  basic  approach  is  to  map  the  exterior  of  an  airfoil  to  the  interior  of 
a  circle.  Polar  coordinates  in  the  mapped  plane  will  be  used.  Spectral 
decomposition  of  the  solution  can  be  used  in  the  mapped  plane.  To  serve  this 
purpose,  the  Euler  equations  in  the  mapped  plane  are  written  in  fully 
conservative  form  by  using  both  physical  and  contravariant  velocities.  These 


equations  are 


3  .  3  .4  3  4 

3t  ^  +  3X  (F)  +  3Y  (G)  =  0 
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(x,y)  are  Cartesian  coordinates  in  the  physical  plane,  (X,Y)  are  polar  coor¬ 
dinates  in  the  mapped  plane,  p  is  the  density,  (u,v)  are  velocity  components, 
(U,V)  are  unsealed  contravariant  velocity  components,  E  is  the  energy,  H  is 
the  specific  enthalpy,  P  is  the  pressure,  J  is  the  Jacobian  of  the  transforma¬ 
tion,  is  the  free-stream  Mach  number,  and  y  is  the  specific  heat  ratio.  P, 
p,  and  the  velocity  vector  (u,v)  are  nondimensionalized  by  their  values  at  the 
free-stream  condition,  E  and  H  are  nondimensionalized  by  the  free-stream  inter¬ 
nal  energy  p  C  T  ,  and  other  variables  at  th?  free-stream  condition  are  com- 

00  y  00 

puted  from  these  variables. 

The  physical  boundary  conditions  are  defined  by  the  solid-wall  condition 
on  the  airfoil  surface  and  the  fact  that  the  disturbances  generated  by  the 
airfoil  propagate  outward  to  infinity.  The  numerical  implementation  of  these 
physical  boundary  conditions  will  be  discussed  later. 


2.2  Numerical  Scheme 

A  computational  mesh  is  created  by  equally  dividing  the  (X,Y)  coordinates 
in  the  mapped  plane.  In  the  mapped  plane,  the  spatial  derivatives  in  X  at 


TN-215/09-84 


r_  •  , 


. 

1 


•  ,s  /  A  A  ••  , 


•' . 


each  mesh  point  are  evaluated  by  application  of  a  fast  Fourier  transform.  The 
derivatives  in  Y  are  evaluated  by  second-order  central  finite  differences. 
Evaluation  of  the  elements  in  the  transformation  matrices  (x^.x^.y^.yy)  is 
performed  in  the  same  manner.  The  singularity  of  the  transformation  at  the 
trailing  edge  is  avoided  by  placing  it  between  two  mesh  points. 

By  using  this  method  of  evaluating  spatial  derivatives,  the  governing 
equations  are  converted  to  a  system  of  ordinary  differential  equations  in  time. 
These  equations  can  be  solved  numerically  by  using  any  of  a  variety  of  well- 
developed  techniques  for  the  solution  of  ordinary  differential  equations. 

An  approximate  fourth-order  Runge-Kutta  scheme  is  used  in  this  work.  The 
algorithm  is  given  by  the  following  equations: 


-►(n)  -*-o  1  -i(n-l) 

q  =  q  +  r^R 


n  =  1, 


k 


(2) 


where  q  represents  the  flow  variables  at  the  mesh  points,  R  represents  the 
terms  with  spatial  derivatives  in  the  equations,  and  n  denotes  the  Runge-Kutta 
step.  For  a  linear  wave  equation,  this  scheme  has  been  shown  to  be  stable  for 
a  CFL  number  less  than  2.8  (Jameson,  Schmidt,  and  Turkel,  1981)  for  a  finite 
difference  scheme,  and  is  stable  for  our  hybrid  scheme  with  a  CFL  number  less 
than  or  equal  to  2.  Following  Jameson  et  al.,  a  local  time  step  that  is  re¬ 
stricted  by  the  CFL  number  is  used.  Because  of  this,  no  physical  interpreta¬ 
tion  should  be  given  to  the  transient  solutions. 


2.3  Filtering 

Filtering  is  required  in  suppressing  the  Gibbs  error.  A  Schumann  filter 
used  by  Gottlieb  et  al.  (1982)  and  given  by  the  following  formula  has  been 
applied  every  35  time  steps  at  a  CFL  number  of  3.5  (see  Section  2.5,  Convergence 
Acceleration) . 

qK  “  0,25(qK+l  +  2qK  +  qK-l);  K  “  (3) 

where  i  and  j  denote  indices  of  the  discrete  points  in  the  X  and  Y  directions, 
respectively.  At  the  shock  wave,  one-sided  filtering  in  the  X  direction  is 
applied  to  preserve  the  sharpness  of  the  shock  wave.  The  Schumann  filter  is 
equivalent  to  a  first-order  artificial  viscosity.  However,  the  filter  is 
applied  only  every  35  time  steps.  The  order  of  the  error  is  higher  than  first 
order.  No  high-mode  smoothing  (Gottlieb  et  al.,  1982)  has  been  applied. 
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We  have  experimented  with  other  filtering  methods,  such  as  derivative 
smoothing  and  artificial  viscosity.  The  former  is  a  low-pass  filter  in  the 
wave  number  space  when  evaluating  derivatives  spectrally.  It  is  expected  to 
yield  a  smoother  residue  and  thus  stabilize  the  calculations.  Numerical 
experiments  do  not  favor  this  approach,  as  the  filter  fails  to  stabilize  the 
computations.  As  to  the  artificial  viscosity,  the  shock  wave  is  smeared  in 
the  computations,  which  defeats  the  nondissipative  nature  of  the  spectral 
scheme. 


2.4  Boundary  Conditions 

The  numerical  implementation  of  boundary  conditions  for  a  hyperbolic 
system  of  partial  differential  equations  is  an  active  research  subject  in 
itself.  Essentially,  on  a  boundary  Y  ,  there  are  four  characteristics  that 
correspond  to  the  speeds  qn,  qt,  qQ  -  c  and  qQ  +  c.  The  respective  charac¬ 
teristic  variables  are  p-c^p,  q^,  p-pcq  ,  and  p  +  p  c  q  ,  where  the 
subscript  o  stands  for  the  quantities  at  the  previous  time  step.  Only  the 
characteristic  variables  carried  on  the  outgoing  characteristics  from  the 
interior  of  the  fluid  domain  can  be  computed  from  the  governing  equations . 

The  characteristic  variables  carried  on  the  incoming  characteristics  must  be 
replaced  by  the  appropriate  boundary  conditions.  The  flow  quantities  can  then 
be  recovered  from  the  combination  of  the  boundary  conditions  and  the  outgoing 
characteristic  variables. 

On  the  solid  surface  Y  ■  0,  there  is  only  one  incoming  characteristic. 

Let  (M,N)  be  the  momentum  along  the  surface  and  normal  to  the  surface, 
respectively,  and  AQ  be  the  symbol  for  the  temporal  change  of  physical 
quantities  Q  as  computed  by  the  Interior  formula,  with  the  subscript  c 
denoting  the  quantities  computed  by  the  interior  formula.  The  following 
formulae  for  the  physical  quantities  at  the  boundary  points  can  then  be  given. 


AP 


AE  +  Y(Y-1)M^  (-  |  AM  +  j  Ap 


P  -  P  +  AP 
c  o 


(4„ 


(5) 


M  /p  +  (AM  -  Ap*M  /p  )/p^ 
o  0  o  o  o 


(6] 
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(7) 
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(8) 
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j 
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Y(y-1)  M2  [M(n)]  /p(n) 

(12) 

.(n) 


where  the  superscript  n  stands  for  the  newly  advanced  quantities,  and  all 
velocities  are  nondimensionalized  by  the  free-stream  velocity. 

At  the  far  field  boundary,  the  treatment  is  essentially  the  same  as  that 
used  by  Jameson  et  al.  (1981)  except  that  the  "extrapolated"  quantities  as 
defined  in  that  paper  are  those  computed  by  the  interior  computations. 

2. 5  Convergence  Acceleration 

To  increase  the  stability  of  the  time-stepping  scheme,  an  additional 
"residue-smoothing"  process  (Jameson  and  Baker,  1983)  has  been  implemented. 

A 'ter  the  residue  R  has  been  evaluated  at  every  mesh  point,  the  residues  are 
smoothed  by  a  linear  transformation  defined  as  follows: 

R  *  (1  -  eS2)-1  (1  ”  g5y)"1  R  (13) 

where  6^  and  6^  are  conventional  finite  difference  operators  in  X  and  Y,  and  e 
is  the  parameter  for  the  residue-averaging  process.  The  new  modified  residue 
field  R  is  used  to  advance  the  solution  in  time.  This  process  alters  the  time- 
dependent  solution  without  changing  its  steady  state.  To  bring  out  the  essen¬ 
tials  of  the  residue-averaging  process,  a  simple  wave  equation  is  considered: 

4>  +  c$  ■  0  .  (14) 

t  x 

The  residue-averaging  process  as  described  is  equivalent,  to  the  lowest 
order,  to  adding  an  additional  term  to  the  original  simple  wave  equation  and 
converting  it  to  the  following  equation: 


A  +  c<J>  -  e(Ax)^  $ 
t  x 


txx 


0 


(15) 


-V- J 


*  -  J 


N*  *■ 
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The  dispersion  relation  for  this  equation  can  be  given  as 

W  o  o  v.-Lb> 

K  1  +  ekz(Axr 

where  (o  is  the  frequency  and  k  is  the  wave  number.  By  increasing  the  para¬ 
meter  e,  the  wave  speed  for  the  high  wave  number  component  is  substantially 
increased.  This  decrease  in  wave  speed  for  the  dangerous  short  waves  contri¬ 
butes  to  the  substantial  increase  in  the  time  step.  In  fact,  Equation  (15)  is 
the  linearized  form  of  a  model  equation  for  long  dispersive  waves  discussed  by 
Benjamin,  Bona,  and  Mahony  (1972),  who  pointed  out  the  numerical  advantage  of 
this  equation  over  the  Korteweg-de-Vries  equation.  Other  means  of  manipulating 
the  dispersion  relation  to  gain  stability  have  been  suggested  (e.g. ,  Gottlieb 
and  Turkel,  1980).  However,  these  methods  do  not  recover  the  original  equation 
in  steady  state,  although  the  error  is  of  higher  order.  The  residue-averaging 
process  substantially  extends  the  stability  boundary.  A  CFL  number  of  3.5  has 
been  used  without  any  difficulty. 
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j.  computed  nesuiis  ioi  nypria  ocneme 

The  initial  effort  employed  a  hybrid  spectral/finite  difference  discretiza¬ 
tion  to  gain  some  experience  with  the  problem.  In  the  circumferential  direc¬ 
tion,  the  variables  are  expanded  in  a  Fourier  series  because  of  the  periodic 
nature  of  the  problem.  In  the  radial  direction,  a  central  finite  difference 
scheme  is  used.  This  hybrid  scheme  is  designed  to  answer  some  questions  in 
applying  the  pseudospectral  method  to  transonic  flows.  The  application  of  the 
pseudospectral  method  to  a  realistic  airfoil  shape  can  be  demonstrated  by  this 
scheme.  Since  we  expect  that  the  shock  wave  will  be  normal  to  the  airfoil 
surface,  the  discontinuity  is  mainly  in  the  X  direction.  Adequate  resolution 
of  the  shock  wave  can  be  achieved  and  the  question  of  convergence  of  the 
spectral  series  can  be  answered  using  the  hybrid  scheme.  Other  properties, 
such  as  the  time-stepping  scheme,  convergence  acceleration,  and  the  filtering 
technique,  can  also  be  studied  with  the  hybrid  scheme.  We  shall  use  a  Karman- 
Trefftz  airfoil  for  this  work  because  of  its  simple  analytical  mapping  from 
the  physical  plane  to  the  interior  of  the  circle.  The  method  can  easily  be 
extended  to  an  airfoil  of  arbitrary  shape  by  using  a  truncated  complex  series 
to  map  the  profile  to  a  circle. 

For  testing  the  solution  algorithm,  flows  around  a  circular  cylinder  are 
computed.  The  pressure  distribution  for  a  subcritlcal  flow  with  *  0.39 
is  given  in  Figure  1.  The  computation  1b  performed  on  a  64x24  grid  (64  points 
in  the  circumferential  direction,  24  points  radially).  It  has  been  computed 
without  filtering.  A  supercritical  case  with  ■  0.45  is  also  computed, 
and  the  results  are  shown  in  Figure  2.  Filtering  is  performed  every  35  time 
steps  with  a  CFL  number  of  3.5.  The  results  agree  with  a  finite  volume  cal¬ 
culation  by  Jameson  et  al.  (1981).  The  shock  wave  has  no  internal  structure 
and  is  sharply  defined. 

A  Karman-Trefftz  airfoil  with  the  following  transformation  from  the  mapped 
plane  £  to  the  physical  plane  z  is  chosen  for  calculations. 


z_  _  (!+!£)*  +  (1-LQ* 
<L  (1+14)*  +  (1-LC)K 


1.9 


L  -  (1  -  rh1/z  -  €  :  C  -  (£  ,n  )  -  (-0.1,0) 

O  O  O  O  0 


d; 

(II 
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Supercritical  non- lifting  flows  with  MWB  0.75  are  computed  on  a 
64x24  grid.  The  results  are  shown  in  Figure  3,  together  with  the  results  from 
a  finite  volume  calculation.  The  hybrid  calculation  shows  again  a  sharply 
defined  shock  wave.  The  agreement  between  the  two  calculations  is  very  good. 
In  particular,  the  positions  of  the  shock  as  defined  by  the  midpoint  of  the 
structure  show  close  agreement.  There  are  discrepancies  Immediately  behind 
the  shock  wave,  however,  and  the  source  of  these  discrepancies  is  not  clear. 
The  pressure  ratio  across  the  shock  wave  using  the  pseudospectral  calculation 
has  been  checked  against  that  using  the  Rankine-Hugoniot  relation  based  on  the 
upstream  Mach  number.  The  error  is  less  than  4  percent. 

To  demonstrate  the  convergence  of  the  Fourier  series,  the  same  case  is 
computed  on  a  32x24  grid.  The  results  are  shown  in  Figure  4,  together  with 
the  results  of  calculations  on  a  denser  mesh  using  a  finite  volume  calcula¬ 
tion.  The  accuracy  of  the  32x24  calculation  is  quite  good.  The  shock 
resolution  of  the  sparse  mesh  calculation  is  comparable  to  that  of  the  64x24 
finite  volume  calculation.  As  expected,  the  finite  volume  calculation  on  the 
sparse  grid  does  not  produce  acceptable  results. 
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4.  Full  Spectral  Scheme 

To  carry  the  work  further,  we  have  attempted  to  construct  a  full  spectral 
method.  Instead  of  using  a  finite  difference  method  for  the  radial  direction 
in  the  transformed  circular  plane,  we  have  used  a  Chebyshev  polynomial 
expansion  for  that  direction. 

We  immediately  encountered  difficulties,  however,  in  using  the  Chebyshev 
spectral  method.  The  first  difficulty  is  the  stability  problem.  Chebyshev 
collocation  points  are  defined  as 


y 


COS  Tt 


CJ-1) 

N 


(19) 


where  N+l  >  J  >  1  is  the  computational  coordinate  and  1  >  y  >  0  is  the  physical 
space.  It  is  easy  to  verify  that  the  spacing  of  the  grid  points  at  the  end 
points  is  asymptotically 


Ay  -  0 


(20) 


For  viscous  problems,  the  end  points  are  usually  on  the  solid  surface,  where 
the  flow  velocity  vanishes.  The  time  steps  can  maintain  a  reasonable  size 
even  though  the  grid  spacing  is  small,  since 


In  the  present  case, 


(21) 

(22) 


where  c  is  the  speed  of  sound.  The  small  grid  size  near  the  end  points  forces 


At  ~ 


(23) 


The  solution  will  take  an  excessively  large  number  of  time  steps  to  develop. 
Also,  because  of  large  variations  in  the  grid  spacing,  the  local  time  step 
approach  used  successfully  in  the  hybrid  method  does  not  seem  to  apply. 

An  attempt  to  stretch  the  Chebyshev  grid  to  achieve  a  more  uniform  grid 
spacing  also  failed.  The  accuracy  near  the  end  points  deteriorates,  which 
affects  the  accurate  application  of  boundary  conditions.  This  deterioration 
of  accuracy  in  a  stretched  grid  can  be  understood  by  taking  the  extreme 
example  of  restoration  to  a  uniform  grid.  Let 

H  ■  j  (1  -  cos  fly)  (24) 
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where  the  computational  space  f|  is  divided  into  Chebyshev  collocation  points 
and  the  corresponding  y  coordinate  is  uniform.  By  the  chain  rule,  we  have 


9__  =  dn  9_ 

3y  dy  3ri 


dn 

dy 


IT 

j  sin  *y 


( 


Since  dri/dy  approaches  zero  at  the  end  points,  accurate  evaluation  of  the 
derivative  at  the  end  points  is  not  possible. 

We  currently  plan  to  investigate  two  other  methods.  The  first  is  to 
discard  the  Chebyshev  method  and  attempt  the  Fourier  polynomial  subtraction 
method  (Gottlieb  and  Orszag,  1977).  The  second  is  to  develop  an  implicit 
method  to  circumvent  the  stability  problem. 
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5 •  Conclusions 

The  hybrid  method  shows  promise  for  the  use  of  spectral  methods  to  accur¬ 
ately  resolve  the  shock  wave  in  a  realistic  problem.  However,  the  filtering 
of  the  solution  certainly  deteriorates  the  accuracy.  Furthermore,  for  complex 
three-dimensional  problems,  it  is  difficult  to  identify  a  supersonic-supersonic 
shock.  Present  means  of  identifying  shock  waves  by  the  sonic  condition  do  not 
apply  there.  Unlike  the  finite  difference  method,  the  residue  of  a  spectral 
method  does  not  decrease  with  the  number  of  time  steps.  At  present,  the 
number  of  supersonic  points  is  used  as  the  indicator  of  convergence. 

The  full  spectral  method  using  a  Chebyshev  expansion  also  presents  diffi¬ 
culties.  Further  investigations  of  alternative  methods  are  required  to  assess 
the  merit  of  a  full  spectral  scheme. 
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Figure  1.  Pressure  Distribution  for  s  Subcritical  Flow  Over  e  Cylinder 
et  Mach  0.39 
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Figure  4.  Comparison  of  Results  with  a  Finite  Volume  Scheme  for 
Karman-Treffte  Airfoil  at  Mach  0.75 
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SUMMARY 

A  hybrid  pseudospectral-f inite  difference  scheme  is  used  to 
calculate  transonic  flow  over  a  two-dimensional  object  using  the 
Euler  equations.  The  exterior  of  the  object  is  mapped  to  the 
interior  of  a  circle.  The  flow  field  variables  are  discre¬ 
tized  using  a  Fourier  series  in  the  circumferential  direction, 
while  a  central  finite  difference  scheme  is  used  in  the  radial 
direction.  We  used  a  four-stage  Runge-Kutta  scheme  including  a 
filter  and  a  residue-smoothing  process.  Transonic  flows  over  a 
circular  cylinder  as  well  as  a  Karman-Tref f tz  airfoil  were 
computed.  The  results  are  compared  to  those  from  finite  volume 
calculations.  It  is  found  that  the  pseudospectral  calculations 
are  able  to  produce  shocks  with  no  internal  structure,  and  fewer 
grid  points  are  needed  to  obtain  the  required,  accuracy. 

INTRODUCTION 

In  recent  years,  there  has  been  strong  interest  in  computa¬ 
tions  of  transonic  flows  using  the  time-dependent  Euler 
equations.  This  interest  stems  in  part  from  the  possibility  of 
shock -generated  vorticity  in  the  flow  field  and  in  part  from  the 
interest  in  numerical  methods  for  a  nonlinear  hyperbolic 
system.  Most  of  the  numerical  methods  are  finite  difference  in 
nature  and  are  second  order  in  accuracy.  To  stabilize  the 
computation  and  to  smooth  the  dispersive  error  for  unsteady 
computations,  either  artificial  dissipative  terms  are  added  to 
the  equations  or  a  built-in  dissipative  mechanism  is  included  in 
the  numerical  scheme.  These  dissipative  terms  usually  cause  the 
shock  wave  to  span  across  three  to  four  grid  points.  To  capture 
a  shock  with  reasonable  accuracy,  one  is  forced  to  use  fairly 
dense  grid  distributions  over  the  region  where  the  shock  wave  is 
expected  to  be. 

The  pseudospectral  method  is  an  alternative  to  the  finite 
difference  method.  It  has  been  applied  successfully  to  many 
smoothly  varying  flows.  The  numerical  analysis  of  the  method 
has  been  given  in  detail  by  Gottlieb  and  Orszag  in  a  monograph 
[1].  Because  of  its  high  rate  of  convergence,  the  method 
usually  requires  relatively  few  terms  of  the  basis  functions  for 
accurate  computations.  In  addition  to  the  spatial  accuracy,  the 
dispersive  error  for  unsteady  computations  is  also  minimized. 

Recently,  efforts  have  been  made  to  apply  the  pseudo¬ 
spectral  method  to  flows  with  shock  waves.  Gottlieb,  Lustman, 
and  Orszag  [2]  have  demonstrated  the  feasibility  of  the  pseudo¬ 
spectral  method  through  the  solution  of  a  one-dimensional  shock 
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tube  problem.  By  using  the  shock-capturing  technique,  they 
showed  that  the  shock  wave  can  be  resolved  within  one  grid  point. 
The  Gibbs  phenomenon  error  due  to  the  discontinuity  can  be  fil¬ 
tered  to  improve  the  accuracy.  Gottlieb,  Lustman,  and  Streett 
[3]  have  attempted  the  two-dimensional  problem  of  the  reflection 
of  an  oblique  shock  from  a  wall.  The  results  from  this  investi¬ 
gation  are  encouraging.  Using  a  fairly  sparse  grid,  they  showed 
that  the  shock  wave  can  be  resolved  within  one  grid  point.  The 
accuracy  of  the  solution  as  compared  to  the  exact  solution  is 
reasonable  considering  the  sparseness  of  the  grid  points. 

In  the  present  work,  we  shall  consider  steady  transonic  flows 
around  an  airfoil  by  solving  the  Euler  equations.  We  shall  ad¬ 
dress  problems  of  applying  the  pseudospectral  method  to  flows 
around  a  complex  geometry,  including  the  development  of  a  time¬ 
stepping  scheme,  and  enhancement  of  the  stability  by  residue 
averaging.  Numerical  experimentation  has  been  used  to  confirm 
convergence  with  a  small  number  of  basis  functions,  and  also  the 
capability  to  treat  shock  waves  with  the  aid  of  filtering. 


GOVERNING  EQUATIONS  AND  BASIC  APPROACH 

The  basic  approach  is  to  map  the  exterior  of  an  airfoil  to 
the  interior  of  a  circle.  Polar  coordinates  in  the  mapped  plane 
will  be  used.  Spectral  decomposition  of  the  solution  can  be 
used  in  the  mapped  plane.  To  serve  this  purpose,  the  Euler 
equations  in  the  mapped  plane  are  written  in  fully  conservative 
form  by  using  both  physical  and  contravariant  velocities.  These 
equations  are 

Jt  +  3X  ^  +  3Y  (^>  “  0  ^ 


V  '  -y£  v  '  H  -[<*  -  DP  *  *]/P  > 

(x,y)  are  Cartesian  coordinates  in  the  physical  plane,  (X,Y)  are 
polar  coordinates  in  the  mapped  plane,  p  is  the  density,  (u,v) 
are  velocity  components,  (U,V)  are  unsealed  contravariant 
velocity  components,  E  is  the  energy,  H  is  the  specific 
enthalpy,  P  is  the  pressure,  J  is  the  Jacobian  of  the  trans¬ 
formation,  Moo  is  the  free-stream  Mach  number,  and  y  is  the 
specific  heat  ratio.  P,  p  and  the  velocity  vector  (u,v)  are 
nondimensionalized  by  their  values  at  the  free-stream  condition, 

E  and  H  are  nondimensionalized  by  the  free-stream  internal  energy 
P«oCvToo,  and  other  variables  at  the  free-stream  condition  are 
computed  from  these  variables. 

The  physical  boundary  conditions  are  defined  by  the  solid- 
wall  condition  on  the  airfoil  surface  and  the  fact  that  the 


disturbances  generated  by  the  airfoil  propagate  outward  to 
infinity.  The  numerical  implementation  of  these  physical 
boundary  conditions  will  be  discussed  later. 

This  initial  effort  employed  a  hybrid  spectral-finite 
difference  discretization  to  gain  some  experience  with  the 
problem.  In  the  circumferential  direction,  the  variables  are 
expanded  in  a  Fourier  series  because  of  the  periodic  nature  of 
the  problem.  In  the  radial  direction,  a  central  finite 
difference  scheme  is  used.  This  hybrid  scheme  is  designed  to 
answer  some  questions  in  applying  the  pseudospectral  method  to 
transonic  flows.  The  application  of  the  pseudospectral  method 
to  a  realistic  airfoil  shape  can  be  demonstrated  by  this  scheme. 
Since  we  expect  that  the  shock  wave  will  be  normal  to  the  air¬ 
foil  surface,  the  discontinuity  is  mainly  in  the  X  direction. 
Adequate  resolution  of  the  shock  wave  can  be  achieved  and  the 
question  of  convergence  of  the  spectral  series  can  be  answered 
using  the  hybrid  scheme.  Other  properties,  such  as  the  time¬ 
stepping  scheme,  convergence  acceleration,  and  the  filtering 
technique  can  also  be  studied  with  the  hybrid  scheme.  We  shall 
use  a  Karman-Tref f tz  airfoil  for  this  work  because  of  its  simple 
analytical  mapping  from  the  physical  plane  to  the  interior  of 
the  circle.  The  method  can  easily  be  extended  to  an  airfoil  of 
arbitrary  shape  by  using  a  truncated  complex  series  to  map  the 
profile  to  a  circle. 


NUMERICAL  SCHEME 


A  computational  mesh  is  created  by  equally  dividing  the 
(X, Y)  coordinates  in  the  mapped  plane.  In  the  mapped  plane, 
the  spatial  derivatives  in  X  at  each  mesh  point  are  evaluated 
by  application  of  a  fast  Fourier  transform.  The  derivatives  in 
Y  are  evaluated  by  second-order  central  finite  differences. 
Evaluation  of  the  elements  in  the  transformation  matrices 
(xx»xy,yx,yy)  is  performed  in  the  same  manner.  The  singularity 
of  the  transformation  at  the  trailing  edge  is  avoided  by  placing 
it  between  two  mesh  points. 

By  using  this  method  of  evaluating  spatial  derivatives,  the 
governing  equations  are  converted  to  a  system  of  ordinary  dif¬ 
ferential  equations  in  time.  These  equations  can  be  solved 
numerically  by  using  any  cf  a  variety  of  well-developed  techni¬ 
ques  for  the  solution  of  ordinary  differential  equations.  An 
approximate  fourth-order  Runge-Kutta  scheme  is  used  in  this 
work.  The  algorithm  is  given  by  the  following  equations: 
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;  n  =  1 
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where  q  represents  the  flow  variables  at  the  mesh  points,  R 
represents  the  terms  with  spatial  derivatives  in  the  equations, 
and  n  denotes  the  Runge-Kutta  step.  This  scheme  has  been  shown 
to  be  stable  for  a  CFL  number  less  than  2.8  [4]  for  a  finite 
difference  scheme,  and  is  stable  for  our  hybrid  scheme  with  a 
CFL  number  less  than  or  equal  to  2.  Following  Jameson,  Schmidt, 
and  Turkel  [4],  a  local  time  step  that  is  restricted  by  the  CFL 
number  is  used.  Because  of  this,  no  physical  interpretation 
should  be  given  to  the  transient  solutions. 


FILTERING 


Filtering  is  required  in  suppressing  the  Gibbs  error.  A 
Schumann  filter  used  by  Gottlieb,  Lustman,  and  Streett  [3]  and 
given  by  the  following  formula  has  been  applied  every  35  time 
steps  at  a  CFL  number  of  3.5  (see  later  section  on  convergence 
acceleration) . 


qK  =  0.25(qK+1  +  2qR  +  qR-1 ) ?  K  =  i,j  (3) 

where  i  and  j  denote  indices  of  the  discrete  points  in  the  X 
and  Y  directions,  respectively.  At  the  shock  wave,  one-sided 
filtering  in  the  X  direction  is  applied  to  preserve  the  sharp¬ 
ness  of  the  shock  wave.  The  Schumann  filter  is  equivalent  to  a 
first-order  artificial  viscosity.  However,  the  filter  is 
applied  only  every  35  time  steps.  The  order  of  the  error  is 
higher  than  first  order.  No  high-mode  smoothing  [3]  has  been 
applied . 


BOUNDARY  CONDITIONS 

The  numerical  implementation  of  boundary  conditions  for  a 
hyperbolic  system  of  partial  differential  equations  is  an 
active  research  subject  in  itself.  Essentially,  on  a  boundary 


lo» 


there  are  four  characteristics  that  correspond  to  the 


speeds  q^  qt,  qn  “  c  and  qn  +  c.  The  respective  characteristic 
variables  are  p  -  c£p ,  qt,  P  -  Poco<Ip'  and  P  +  Poco<Jn»  where 
the  subscript  o  stands  for  the  quantities  at  the  previous  time 
step.  Only  the  characteristic  variables  carried  on  the  outgoing 
characteristics  from  the  interior  of  the  fluid  domain  can  be 
computed  from  the  governing  equations.  The  characteristic 
variables  carried  on  the  incoming  characteristics  must  be 
replaced  by  the  appropriate  boundary  conditions.  The  flow 
quantities  can  then  be  recovered  from  the  combination  of  the 
boundary  conditions  and  the  outgoing  characteristic  variables. 

On  the  solid  surface  Y  =  0,  there  is  only  one  incoming 
characteristic.  Let  (M,N)  be  the  momentum  along  the  surface 
and  normal  to  the  surface,  respectively,  and  AQ  be  the  symbol 
for  the  temporal  change  of  physical  quantities  Q  as  computed  by 
the  interior  formula,  with  the  subscript  c  denoting  the 
quantities  computed  by  the  interior  formula.  The  following 
formulae  for  the  physical  quantities  at  the  boundary  points  can 
then  be  given. 


AP  =  AE  +  y(y-1)M„ 


l  M  AM  ^  1  M2  *  \ 

p  AM  +  2  ^2  t0 J 


+  AP 
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M0/P0  +  -  Ap.Mo/po)/po 


pc  +  Ap 


P(n)  «  P  -  ym2 
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p(n)  "  pc  +  “^2  (p(n)  "  Pc)/co 

c  ym 2  C  ° 

M(n)  =  p (n) . (M/p)c 

N*n*  *  0 

E(n)  =  p(n)+  1  Y(Y_X)  m2 

where  the  superscript  n  stands  for  the  newly  advanced  quantities, 
and  all  velocities  are  nondimensionalized  by  the  free-stream 
velocity. 

At  the  far  field  boundary,  the  treatment  is  essentially  the 
same  as  that  used  by  Jameson,  Schmidt,  and  Turkel  [4]  except 
that  the  "extrapolated"  quantities  as  defined  in  that  paper  are 
those  computed  by  the  interior  computations. 


CONVERGENCE  ACCELERATION 

To  increase  the  stability  of  the  time-stepping  scheme,  an 
additional  "residue-smoothing"  process  [5]  has  been  implemented. 
After  the  residue  R  has  been  evaluated  at  every  mesh  point,  the 
residues  are  smoothed  by  a  linear  transformation  defined  as 
follows : 


(10) 

(11) 

(12) 


R  *  (1  -  efi2)"1  (1  -  efiy)'1  R 


(13) 


where  fix  and  are  conventional  finite  difference  operators 
in  X  and  Y,  and  e  is  the  parameter  for  the  residue-averaging 
process.  The  new  modified  residue  field  R  is  used  to  advance 
the  solution  in  time.  This  process  alters  the  time-dependent 
solution  without  changing  its  steady  state.  To  bring  out  the 
essentials  of  the  residue-averaging  process,  a  simple  wave 
equation  is  considered: 

$t  +  c$x  *  0  *  (*4) 


The  residue-averaging  process  as  described  is  equivalent, 
to  the  lowest  order,  to  adding  an  additional  term  to  the 
original  simple  wave  equation  and  converting  it  to  the 
following  equation: 

4>t  +  c4>x  -  eUx)2  $txx  *  0  .  (15) 

The  dispersion  relation  for  this  equation  can  be  given  as 


U) 

k 


_ c _ 

1  +  ek2(fix)2 


(16) 


where  w  is  the  frequency  and  k  is  the  wave  number.  By  increasing 
the  parameter  c ,  the  wave  speed  for  the  high  wave  number  com¬ 
ponent  is  substantially  increased.  This  decrease  in  wave  speed 
for  the  dangerous  short  waves  contributes  to  the  substantial 
increase  in  the  time  step.  In  fact,  Equation  (15)  is  the 
linearized  form  of  a  model  equation  for  long  dispersive  waves 


discussed  by  Benjamin,  Bona,  and  Mahony  [6],  who  pointed  out 
the  numerical  advantage  of  this  equation  over  the  Korteweg-de- 
Vries  equation.  Other  means  of  manipulating  the  dispersion 
relation  to  gain  stability  have  been  suggested  (e.g.,  [7]). 
However,  these  methods  do  not  recover  the  original  equation  in 
steady  state,  although  the  error  is  of  higher  order.  The 
residue-averaging  process  substantially  extends  the  stability 
boundary.  A  CFL  number  of  3.5  has  been  used  without  any 
difficulty. 


COMPUTED  RESULTS 

For  testing  the  solution  algorithm,  flows  around  a  circular 
cylinder  are  computed.  The  pressure  distribution  for  a  sub- 
critical  flow  with  Moo  =  0.39  is  given  in  Figure  1.  The  com¬ 
putation  is  performed  on  a  64x24  grid  (64  points  in  the 
circumferential  direction,  24  points  radially).  It  has  been 
computed  without  filtering.  A  supercritical  case  with 
Moo  =  0.45  is  also  computed,  and  the  results  are  shown  in 
Figure  2.  Filtering  is  performed  every  35  time  steps  with  a 
CFL  number  of  3.5.  The  results  agree  with  a  finite  volume 
calculation  by  Jameson,  Schmidt,  and  Turkel,  [4].  The  shock 
wave  has  no  internal  structure  and  is  sharply  defined. 

A  Karman-Tref f tz  airfoil  with  the  following  transformation 
from  the  mapped  plane  £  to  the  physical  plane  z  is  chosen  for 
calculations . 

z  =  tKmlm9  (17) 

Kh  (l+LO*  +  (l-LO* 

L  =  (1  -  n2)1/2  -  *  ^c'T1o)  =  (-0.1,0)  (18) 

Supercritical  non-lifting  flows  with  Moo  =  0.75  are  computed 
on  a  64x24  grid.  The  results  are  shown  in  Figure  3,  together 
with  the  results  from  a  finite  volume  calculation.  The  hybrid 
calculation  shows  again  a  sharply  defined  shock  wave.  The 
agreement  between  the  two  calculations  is  very  good.  In 
particular,  the  positions  of  the  shock  as  defined  by  the 
midpoint  of  the  structure  show  close  agreement.  There  are 
discrepancies  immediately  behind  the  shock  wave,  however.  The 
source  of  these  discrepancies  is  not  clear.  The  pressure  ratio 
across  the  shock  wave  using  the  pseudospectral  calculation  has 
been  checked  against  that  using  the  Rankine-Hugoniot  relation 
based  on  the  upstream  Mach  number.  The  error  is  less  than 
4  percent. 

To  demonstrate  the  convergence  of  the  Fourier  series,  the 
same  case  is  computed  on  a  32x24  grid.  The  results  are  shown 
in  Figure  4,  together  with  the  results  of  calculations  on  a 
denser  mesh  using  a  finite  volume  calculation.  The  accuracy  of 
the  32x24  calculation  is  quite  good.  The  shock  resolution  of 
the  sparse  mesh  calculation  is  comparable  to  that  of  the  64x24 
finite  volume  calculation.  As  expected,  the  finite  volume 
calculation  on  the  sparse  grid  does  not  produce  acceptable 
results . 
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Figure  1.  Pressure  distri¬ 
bution  for  a  subcritical 
flow  over  a  cylinder  at 
Mach  0.39 


"A 


Figure  3.  Pressure  distri¬ 
bution  on  a  Karoan-Tref f tz 
airfoil  at  Mach  0.75  using 
a  64x24  grid  (350  steps) 


Figure  2.  Pressure  distri> 
bution  for  a  supercritical 
flow  over  a  cylinder  at 
Mach  0.45 
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Figure  4.  Comparison  of 
results  with  a  finite 
volume  scheme  for  a  Karman- 
Trefftz  airfoil  at  Mach  0.75 
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CONCLUSIONS  AND  ACKNOWLEDGEMENTS 


Several  conclusions  can  be  drawn  from  the  present 
investigation. 

(1)  In  computing  flows  with  shock  waves,  the  Gibbs  error  can  be 
filtered  to  produce  accurate  results.  Because  of  the  rapid 
convergence  of  the  Fourier  series,  fewer  grid  points  are 
required  than  with  the  lower  order  difference-type  scheme. 

(2)  A  shock  wave  without  internal  structure  can  be  produced. 
This  capability  also  contributes  to  the  accuracy  of  the 
method  in  that  fewer  grid  points  are  required  to  resolve 
the  shock  wave. 

(3)  Application  of  the  pseudospectral  method  to  flows  around  a 
realistic  geometry  is  possible  using  a  mapping  technique. 

The  authors  are  indebted  to  Mr.  Morton  Cooper  of  Flow 
Research  Company  for  suggesting  the  problem.  Discussions  with 
Professor  Steven  Orszag  have  also  been  very  helpful.  This  work 
was  supported  by  the  Air  Force  Office  of  Scientific  Research 
and  the  Office  of  Naval  Research  under  Contract  No. 
F49620-82-C-0022 . 
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1 .  Introduction 

Work  is  in  progress  at  Flow  Industries  on  the  direct  numerical  simulation 
of  complex  VTOL  flows  using  the  full  three-dimensional,  time-dependent 
Navier-Stokes  equations.  The  objective  of  this  numerical  simulation  is  to 
compute  accurately  the  details  of  the  flow  field  and  to  achieve  a  better 
understanding  of  the  physics  of  the  flow,  including  the  role  of  initial 
turbulence  in  the  jet,  the  influence  of  forward  motion  on  hover  aerodynamics, 
the  collision  zone  and  fountain  characteristics,  and  the  jet  structure  and 
entrainment  process  in  the  transitional  flight  regime.  The  results  of  this 
work  can  be  used  to  evaluate  the  merit  of  various  models  suggested  in  the  past 
or  can  be  used  to  construct  a  new  model.  This  work  will  also  allow  the 
assessment  of  wing-jet-ground  interference  effects  and  the  accurate  prediction 
of  their  associated  forces  and  moments,  which  is  required  for  the  design  and 
optimization  of  VTOL  aircraft.  This  note  describes  the  work  completed  at  Flow 
in  the  second  year  of  a  program  in  which  VTOL  aerodynamics  are  being 
investigated  numerically. 

The  problem  under  investigation  is  that  of  an  infinite  row  of  jets 
impinging  on  the  ground  (see  Figure  1).  This  problem,  which  contains  the 
essential  features  of  twin  jets  impinging  on  the  ground  (see  Figure  2), 
simulates  the  hovering  configuration.  The  choice  of  a  row  of  jets  provides 
the  periodic  property  of  the  flow  field,  which  allows  approximation  of  the 
flow  properties  in  the  periodic  direction  by  a  truncated  Fourier  series.  The 
spectral  method  may  therefore  be  used  in  the  periodic  direction,  while  finite 
difference  approximations  are  used  in  the  vertical  z  direction  and  the 
y  direction  normal  to  the  row  of  jets.  The  jets  may  be  inclined  in  the 
y  direction,  which  leads  to  a  configuration  associated  with  an  aircraft  in 
pitch  while  hovering.  By  imposing  a  cross  flow  in  the  y  direction,  it  is 
possible  to  study  the  effects  of  the  aircraft's  forward  motion  during  takeoff 
and  transition. 

A  computer  code  that  solves  the  time-dependent  Navier-Stokes  equations  has 
been  developed  with  the  purpose  of  numerically  simulating  the  problem  of  an 
infinite  row  of  jets  impinging  on  the  ground.  The  code  presently  uses  finite 
difference  approximations  in  all  three  spatial  directions,  and  it  uses  a 
first-order  time-differencing  scheme.  Modifications  are  in  progress  that  will 
allow  the  use  of  the  spectral  method  in  the  periodic  direction  and  the  use  of 
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a  second-order  time-differencing  scheme.  Subgrid-scale  modeling,  which  allows 
the  solution  of  problems  at  high  Reynolds  numbers,  will  also  be  introduced 
into  the  code.  Although  the  code  is  not  in  its  final  form,  it  has  been  used 
to  obtain  solutions  that  indicate  the  main  features  of  VTOL  aerodynamics. 

In  this  note,  the  governing  equations  and  the  boundary  conditions  used  in 
the  code  are  summarized  in  Section  2.  The  method  of  solution  is  discussed  in 
Section  3,  and  preliminary  examples  of  solutions  using  the  code  are  presented 
in  Section  4. 


m 
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2.  Governing  Equations 

The  governing  equations  are  the  Navier-Stokes  equation 

<?t  +  (q*?)q  ■ _7p  +  ib; 72 2,  (1 

and  the  continuity  equation 

7*q  -  0  (2 

where  q  is  the  velocity  vector  and  p  is  the  pressure.  When  the  Reynolds 
number  is  too  large  to  resolve  numerically  the  entire  range  of  energetic  scales 
filtering  will  be  used  to  eliminate  the  smaller  (subgrid-scale)  motions. 
Filtering  Equation  (1)  introduces  new  terms,  similar  to  Reynolds  stress  terms 
obtained  in  the  Reynolds-averaged  equations,  that  contain  the  effect  of  the 
subgrid-scale  motions  on  the  numerically  resolved  motions.  We  plan  to  use 
standard  procedures  to  handle  these  terms  and  will  introduce  them  into  our 
numerical  scheme  at  a  later  date.  By  taking  the  divergence  of  Equation  (1), 
the  following  Poisson  equation  governing  the  pressure  is- obtained: 

V2p  -  -7*[(q*V)q]  -  (V‘q)t  +  ^  V2(V«q)  .  (3 

Substituting  Equation  (2)  into  Equation  (3)  leads  to  the  following  pressure 
equation: 

V2p  -  -V«[(q*V)q]  .  (4 


The  system  of  Equations  (1)  and  (4)  is  equivalent  to  the  original  system. 
Equations  (1)  and  (2),  and  is  used  here  instead  of  the  original  set  of 
equations.  The  vector  equation  (1)  is  solved  subject  to  the  periodicity 
condition  in  the  x  direction;  a  weak  outflow  condition 


3u  _  3w 
9y  ciy 


0 


» 


rfith  3v/3y  being  determined  from  the  continuity  equation,  is  applied  at 
the  side  boundaries  of  the  computational  domain  (y  ■  yfi,  y  ■  y^),  and  a 
no-slip  condition  is  applied  at  the  bottom  boundary,  z  ■  z  ,  and  the  top 


(5 
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boundary,  z  ■  z  ,  outside  the  jet  region.  In  the  jet  region,  the  inflow 
condition 

q(x,y,za)  -  f(x,y) 

is  specified.  Equation  (4)  is  solved  subject  to  the  condition 


at  the  side  boundaries  and  subject  to  the  Neumann  boundary  conditions 
determined  from  Equation  (1)  at  the  top  and  bottom  boundaries. 
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3.  Method  of  Solution 


The  finite  difference  approximations  to  Equations  (1)  and  (A)  are  written 
at  the  mesh  points  of  a  staggered  grid.  For  advancing  the  solution  from  time 
tn  to  time  tn+\  where  tn  *  nAt,  tn+^  =  (n+1)  At  and  At  is  the  time  step,  the 
following  first-order  scheme  is  used: 

qn+1  -  qU  +  At[-(qn*7)qn  -  VpD  +  ^  V2qn]  (8! 


where  pQ  is  determined  so  that  mass  conservation  is  assured  at  time  tn 

It  is  determined  by  solving  the  finite  difference  approximation  to  the  equation 

V2pn  =  -V*[(qn-7)qn]  .  (9 

This  equation  is  solved  by  using  a  direct  (noniterative)  fast  Laplace  equation 
solver. 
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4 .  Numerical  Examples 

The  examples  presented  here  are  preliminary  examples  that  have  been  solved 
using  the  developed  computer  code.  A  relatively  coarse  numerical  mesh  is  used 
and  the  Reynolds  number  is  assumed  to  be  low  enough  so  that  filtering  is  not 
required.  The  results  presented  here  are  not  intended  to  be  an  accurate 
simulation  of  VTOL  flow  configurations.  Nevertheless,  they  do  indicate  the 
main  features  of  these  flows. 

For  all  the  examples  presented  here,  the  plane  y  =  y^  is  assumed  to  be  a 
plane  of  symmetry.  Unless  otherwise  stated,  the  computational  domain  is 
defined  by  (see  Figure  1) 

0  =  x  .  <x  <  x ,  =  1 
j  -  -  f 

-2  -  yB  <  y  iyb  =  2 

0  =  Zg  <  z  <  za  =  1 


where  all  dimensions  are  normalized  by  the  jet  diameter.  The  jet  velocity 
profile  in  the  direction  of  the  jet  axis  is  assumed  to  be  given  by 


Qj(r) 


where  Rj  is  the  jet  radius,  r  is  the  distance  from  the  jet  axis  and  velocities 
are  normalized  by  the  jet  maximum  velocity.  The  Reynolds  number  in  the 
examples  is  based  on  the  jet  diameter  and  the  maximum  jet  velocity. 


Example  1: 

In  this  example  the  jet  axis  is  assumed  to  be  normal  to  the  ground  plane 
(a  *  90°)  and  there  is  assumed  to  be  no  cross  flow  (V  =  0).  The  Reynolds 
number  is  given  by  Re  =  300.  An  18x72x18  (x,y,z)  mesh  is  used. 

Figures  3  through  12  show  the  main  features  of  the  flow  generated  by  a  row 
of  vertical  jets  impinging  on  the  ground.  The  velocity  vectors  in  the  planes 
x  =  Xj,  x  =  Xj  are  shown  in  Figures  3  and  4,  respectively.  The  fan-shaped 
fountain  that  results  from  the  collision  of  the  two  wall  jets  is  apparent  in 
Figure  4.  The  jet,  the  wall  jet  and  the  fountain  are  apparent  in  Figure  5. 

In  both  Figures  5  and  6  a  downward  motion  exists  as  the  plane  x  =  x^  is 
approached,  while  an  upward  motion  exists  as  the  plane  x  =  x^  is  approached. 
However,  the  relative  magnitudes  of  these  motions  are  reversed  in  the  two 
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figures.  Figure  5  indicates  that  the  downward  motion  in  the  jet  is  greater 
than  the  upward  motion  in  the  fountain.  The  plane  of  Figure  6  does  not  pass 
through  the  jet.  There,  the  fountain  upward  motion  is  relatively  stronger 
than  the  downward  motion  at  x  ■  x^.  Figures  7  through  10  are  pressure 
contours  that  indicate  high-pressure  areas  in  the  zones  of  jet-ground 
impingement,  wall  jet-wall  jet  collision  and  fountain  impingement  on  the  upper 
boundary.  Figures  11  and  12  are  contour  plots  for  the  vorticity  components. 

Example  2; 

In  this  example  the  jet  axis  is  assumed  to  be  inclined  at  an  angle 
a  *  60°  to  the  ground.  A  cross  flow  of  V  =  0.2  is  also  assumed.  The 
Reynolds  number  is  given  by  Re  =  300.  An  18x72x18  (x,y,z)  mesh  is  used. 

Figures  13  through  22  show  the  main  features  of  the  flow  generated  by  a 
row  of  inclined  jets  impinging  on  the  ground  in  a  cross  flow.  In  Figure  13 
the  ground  vortex  formed  by  the  interaction  of  the  cross  flow  and  the  wall  jet 
is  apparent.  The  effect  of  the  cross  flow  on  the  fan-shaped  fountain  is  shown 
in  Figure  14,  where  it  is  no  longer  symmetric. 

For  the  problem  of  a  jet  in  a  cross  flow,  two  basic  configurations  are 
relevant  to  VTOL  aerodynamics.  In  the  first  configuration,  the  jet  impinges 
on  the  ground.  The  main  features  of  this  flow  are  indicated  in  Example  2.  A 
second  configuration  results  as  the  distance  between  the  aircraft  and  the 
ground  becomes  large  and/or  as  the  forward  aircraft  speed  becomes  large.  In 
this  case,  the  jet  does  not  impinge  on  the  ground.  This  configuration  is  used 
in  the  following  example. 

Example  3; 

In  this  example  a  =  90°,  V  =  0.7  and  Re  =  60.  A  7x28x14  mesh  is  used. 

The  computational  domain  is  defined  by 

0=x.<x<x“l 
J  -  ~  f 

-2  -  h  i  y  i  h  m  1 

0  *  z  <z<z  *2 
g  -  -  a 

Figures  23  through  32  show  the  main  features  of  this  flow.  Figure  23 
indicates  that  the  jet  changes  its  direction  before  it  reaches  the  ground.  As 
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indicated  in  Figure  24,  no  fountain  flow  develops  in  this  example  since  there 
are  no  wall  jets.  The  double  vortex  generated  by  the  jet-cross  flow  interac¬ 
tion  is  shown  in  Figure  26.  As  indicated  by  the  pressure  contours  shown  in 
Figure  27,  a  high-pressure  region  develops  upstream  of  the  jet,  while  a 
low-pressure  region  develops  downstream  of  the  jet  in  its  wake. 

Example  4: 

This  example  is  the  same  as  Example  3  with  the  exception  that  V  *  0. 
Figures  33  through  42,  which  show  the  main  features  of  the  flow  for  this  case, 
thus  allow  comparison  between  the  zero  cross-flow  case  (given  here)  and  a 
cross-flow  case  (Example  3)  for  the  same  configuration. 

Example  5: 

This  example  is  similar  to  Example  1  (a  *  90°,  V  *  0).  However,  a 
relatively  fine  computational  mesh  is  used  here,  which  allows  the  use  of  a 
relatively  large  Reynolds  number  (Re  =  600) .  Symmetry  is  assumed  in  the  x 
direction  in  addition  to  the  y  direction.  Therefore,  the  computational  domain 
is  defined  by 

0*Xj_<x£Xj  =  l 

0  -  »j  i  y  £  yb  '  2 

0  =  z  <  z  <  z  =1 
g  -  -  a 

A  24x48x24  computational  mesh  is  used. 

The  results  of  this  calculation  are  indicated  in  Figures  43  through  52. 
Qualitatively  similar  effects  to  those  observed  in  the  first  example  are 
indicated.  However,  certain  effects,  such  as  the  propagation  of  the  initial 
vortex,  while  observed  for  the  high-Reynolds-number ,  fine-mesh  calculation 
(see  Figure  53)  are  not  observed  for  the  low-Reynolds-nuraber,  coarse-mesh 
calculation  (see  Figure  54). 
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Figure  2.  Schematic  of  Fiow  Field  About  a  Hovering  VTOL  Aircra 
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Figure  6.  Example  1:  Velocity  Vectors  in  the  Plane  y  = 
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Figure  7.  Example  1:  Preaeure  Contours  in  the  Plane  x 
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Figure  8.  Example  1:  Pressure  Contours  in  the  Plane 
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Figure  9.  Example  1 :  Pressure  Contours  in  the  Plane  z  =  zg 


Figure  10.  Example  1:  Pressure  Contours  in  the  Plane  z  =  za 
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Figure  13.  Example  2:  Velocity  Vectors  in  the  Plane  x  =  x 
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Figure  14.  Example  2:  Velocity  Vectors  in  the  Plane  x  =  Xf 
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Figure  15.  Example  2:  Velocity  Vectors  in  the  Plane  y  *  yj 
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Figure  16.  Example  2:  Velocity  Vectors  in  the  Plane  y  =  y  ,| 
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Figure  26. 


Example  3:  Velocity  Vectors  In  the  Plane  y  ■ 
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Figure  30.  Example  3:  Pressure  Contours  in  the  Plane  z  =  za 


Figure  33.  Example  4:  Velocity  Vectors  in  the  Plane  x  =  xj 


Figure  34.  Example  4:  Velocity  Vectors  in  the  Plene  x  «  Xf 
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Figure  36.  Example  4:  Velocity  Vectora  in  the  Plane  y  =  yj 


Figure  37.  Example  4:  Pressure  Contours  in  the  Plane 
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Figure  45.  Example  5:  Velocity  Vectors  In  the  Plane  y  ■  yj 
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Figure  49.  Example  5:  Preaeure  Contours  in  the  Plane  z 


Figure  51.  Example  5:  x-Vorticity  Component  in  the  Plane  x  -  xj 
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